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A recently developed continuous time solver based on an expansion in hybridization about an 
exactly solved local limit is reformulated in a manner appropriate for general classes of quantum 
impurity models including spin exchange and pair hopping terms. The utility of the approach 
is demonstrated via applications to the dynamical mean field theory of the Kondo lattice and 
two-orbital models. The algorithm can handle low temperatures and strong couplings without 
encountering a sign problem. 

o 

(N ■ I. INTRODUCTION 

One of the fundamental challenges of theoretical condensed matter physics is the accurate solution of quantum 
impurity models. These, in general terms, consist of a Hamiltonian involving a finite number of states and a hy- 
bridization process which allows particle exchange with one or more "reservoirs" of particles. They are important 
both in their own right and as a crucial ingredient in the dynamical mean field (DMFT) method of approximating 
the properties of interacting fermions on a lattice. Examples include the familiar Kondo and Anderson Hamiltonians 
and their generalization to multi-spin and multi-orbital cases, as well as to the "embedded plaquettes" used in the 
recently developed cluster extensions of dynamical mean field theory HHII. 

Quantum impurity models may be formulated as quantum field theories in zero space and one time dimension, 
\ and the reduced dimensionality suggests that numerical approaches should be feasible. However, up to now general 

■ quantum impurity models have to a large degree resisted numerical attack. A special but conceptually crucial model, 
the one-orbital Anderson impurity model, has been studied in detail but the techniques (the Hirsch-Fye discrete 
Hubbard- Stratonovich transformation [j| and exact diagonalization which work relatively well in this case have 
proven difficult to extend to wider classes of models of physical interest. 

Cn| . One issue is that the Hirsch-Fye method cannot easily be applied to models with interactions other than direct 
density-density couplings. In particular, there is no good decoupling for the exchange and "pair hopping" terms which 
are important to the physics of partially filled d-levels. A scheme proposed by Sakai et al. @ has been used in some 
DMFT studies |ialll|. but the method has a severe sign problem which prevents calculations at low temperatures. 
Another issue with Hirsch-Fye and similar methods is time discretization, and in particular the fine grid spacing 
required to capture the short time behavior of the Green function. The computational burden in Hirsch-Fye type 

\^ " methods grows as the cube of the (large) grid size, which must be increased linearly with interaction strength and 

f^) , inverse temperature. This severely restricts the accessible parameter range. 

■ The exact diagonalization method [|| represents the continuous density of states of the reservoir by a small number 
of levels-but the number of levels required scales linearly with the number of orbitals included while the computational 
burden grows exponentially with the number of levels. This limits the applicability of the method to models with a 
small number of orbitals, although some results have been presented for three orbital models and four-site clusters 
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Recently, a new class of impurity solvers has been developed [13, [Hi, based on the stochastic evaluation of a 
diagrammatic expansion of the partition function. Two complimentary approaches are possible, based on a weak- 
coupling expansion in powers of the coupling constants [T^ or an expansion in powers of the impurity-bath mixing 
|13|. These algorithms, which require neither auxiliary fields nor a time discretization, have been shown to provide 
considerable improvements over the Hirsch-Fye method for the one-orbital Anderson model. The weak coupling 
approach has also been successfully applied to more complicated models ^ij, and an interesting hybrid scheme 
involving a Hirsch-Fye decoupling of density channel interactions and an expansion in exchange interactions has very 
recently been applied to multiorbital models 0, ^| . 

In Ref. ^3 we have demonstrated the usefulness of the hybridization expansion approach for the single site Hubbard 
model. Its power relies on the fact that the order of perturbation which is needed decreases as the interaction strength 
increases. The algorithm was found to allow access to extremely low temperatures, even in the presence of strong 
interactions. But the formulation given in Ref. fl3l ] was specific to models (such as the Hubbard model) with only 
density-density interactions. In this paper we present a matrix formulation which generalizes the method to wide 
classes of impurity models. To demonstrate the power of the hybridization expansion approach we use it to calculate 
physical properties of the dynamical mean field approximation to the Kondo lattice model (for which only very few 
DMFT calculations have been attempted) and the multiorbital Anderson model. 
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II. FORMALISM 



A general impurity model contains fermions labeled by quantum numbers a = 1,...,N (denoting for example 
site, spin and orbital indices), interacting with each other, coupled to local degrees of freedom T (representing for 
example spin or phonon fields) and hybridized with "bath" fermions. The latter have a continuous density of states 
which we parametrize by "momentum" p. It is convenient to assemble the fermion fields and the bath fermions into 
TV-component spinors ip and b, respectively. The general Hamiltonian is then 

H = i?loc + -Hbath + -Hhyb + -^hyb> (•"•) 

with 

H loc = ^Qip-T + h.c. + H T + ]T U abcd ijl^ c tp d , (2) 

-Hbath = ^e p 6 p 6p, (3) 
p 

H hyh = $>V p 6t. (4) 

p 

We have assumed here that the fermion-fermion interaction is of the conventional four-fermion type, but the extension 
to more general forms is immediate. Similarly, we have assumed a bilinear coupling (specified by some matrix Q) 
between the local fermions and the spin and lattice degrees of freedom represented by T, but more general interactions 
are easily included. 

The "bath" fermions are assumed to be orthogonal and to have free fermion correlations while V is an N x N 
hybridization matrix, which has to be determined in a self-consistent manner. In the impurity models known to us it 
is possible to find a representation in which -ffbath an d V are simultaneously diagonal, that is 

tfhyb + H hath = + E = E H W> + E ^bath, (5) 

a.p a,p a a 

and we make this assumption throughout the rest of this paper. 
The impurity model partition function Z may then be expressed as 

Z = Z h ^T H (T T e- ft ^^oo(r)+H bath (r)+E tl (^ yb (r) + < b (r))^ > (g) 

with Z bath = Tne-? 11 ^ and (.) 6 = Tr b [.]/Z hath . 

We expand Eq. © in the hybridizations ipaVpb 1 ^ and b p V p a * Each term in the expansion must have the same 
number of i[> a and operators, so 

Z = ^J^^^^-W+^WnE^), (7) 

a k a 



rP rP rP rP rP rP 

E E K v pT-K a v pV d n / *n- d < ^•••/ < 

oi,...,Pk a p'...,pi J0 Jti J ^ a -i Jo Jt 'i Jt L-i 



Zk a — 

Pl,---,Pk a p' lt ...,p' ka 

^MriKlir^b;,^^^ K>1«J, (8) 



where we have used the l/k a \ in time ordering the tp's and i/j^s. We now take the expectation value over the bath 
states. The unprimed and primed p-indices must always occur in pairs pi ~ p'j = P and tracing over the bath states 

thus yields a factor {V^e-^h^ /(e-^" + 1) if rj > n and \V p a \ 2 e~ t " { - T '^ +p) / {e' 13 ^ + 1) if rj < n . By defining 
the hybridization function F a {r) as 

P ( , _ / E P \V p ^e-^~T)/(e-^ + 1) r > 

~ I E P -\V p a \ 2 e-^-^/(e~^ + l) r < ' W 

the expectation value of the ^-operators can be expressed as the determinant of a matrix M~ x with elements 

M- 1 (i,j) = F a (r i -r;). (10) 
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FIG. 1: (Color online) Every Monte Carlo configuration can be represented by a sequence of operators on the time interval 
< t < /3 (we let time run from right to left to be consistent with the time ordering convention). Different colors correspond 
to different flavors, while full (empty) circles represent creation (annihilation) operators. The Monte Carlo moves consist of 
random insertions or deletions of pairs of operators in the different channels. 



Note that 



F(-iu n ) 



J dTe- l ^ T F{r) = J dujJ2\ V P 



,S(uj - e p ) 



(11) 



so that the hybridization functions F are the same as those defined in Ref. 
defined "Weiss function" Q^ 1 [| by F(~iuj n ) = iuj n + \i — Q^ x (iio^). 
The partition function finally becomes 



[l3j and are related to the conventionally 



Z = 



ZbathTr^ 

-a 



a k a 

rP rP fP [0 rP rP 

/ dn dr 2 ... dr ka / dr{ dr' 2 ... I dr' ka de^M" 1 )*. 

JO J Ti J Tk a — i JO i/r{ T L —l 



(12) 



(13) 



with s a a sign determined by the signature of the permutation which permutes the a-flavored field operators from 
their time-ordered sequence (smallest r shifted to the right) into the alternating order tpa^i)^ (T[)ip a (T 2 )ip a (r^) ■ ■ ., 
and st t compensating for an eventual sign change produced by the time ordering of all the operators. The sign 

e 

factor s a arises from the /3-antiperiodic definition of F a (Eq. ©) and is the generalization of the signs denoted <5 r s 
in Ref. The sign st t is merely a consequence of the notation in Eq. (T^J, where we grouped together all of the 
operators corresponding to a given flavor a. If all the ip operators, irrespective of flavor, are placed in the order in 
which they occur, there is no additional sign. 



III. MONTE CARLO PROCEDURE 



Eqs. H12fl and l|13|l show that the partition function may be expressed as a sum over configurations consisting of 
2n = 2^ a fc a operators {Oi(Ti)}o<Ti<T 2 <...<T 2 „</3- Of these operators, k a are creation operators tp'l and another k a 
are destruction operators ip a and they are connected in all possible ways by hybridization functions F a (this is the 
interpretation of the determinant). Sandwiched in between the O's arc time evolution operators -K\ c, defined as 

K 1oc (t) = e~ H ^\ (14) 

A typical configuration can thus be illustrated by a sequence of dots on an interval [0, (3) representing imaginary 
time (see Fig. Each color corresponds to a different flavor a, while full and empty dots represent creation and 
annihilation operators, respectively. The weight of such a configuration is given by 

w({Oi{n)}) - Tr [K loc (j3 - r 2 „)0 2 „(r 2n ) . . . 2 (t 2 )A- 1oc (t 2 - t x )O x (tx)K 1oc {t x )\ dr x . . . dr 2n JJ (det M~ x )s a . (15) 

a 

A Monte Carlo procedure which samples the whole configuration space is obtained by randomly inserting and 
removing pairs of operators in the a channel (a = 1, . . . , N), or changing their position on the time interval. The 
detailed balance condition for insertion/removal of a pair in channel a reads 

p({Q}2a) p({Q} 2 n+ 2 ) _ P 2 Tr[K loc ((3 - r 2n+2 )6 2ra+2 (f 2 » +2 ) . . . Q 1 (f 1 )A^ loc (f 1 )] detM- 1 ^ 
p({6} 2n+2 ) ^ p({0} 2n ) (fca + 1) 2 Tr[K loc ((3-T 2n )0 2n (T 2n )...0 1 (T 1 )K loc (T 1 )} defM," 1 ^' 

and can be satisfied for example by using the Metropolis algorithm. In each update, it is therefore necessary to 
compute both the determinant of the new i^-matrix, detM" 1 , and the trace of the new sequence of field operators 
and propagators. This latter task is simplified by writing all the operators in the eigenbasis of H\ oc . 
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FIG. 2: (Color online) If the creation and annihilation operators for each flavor must occur in alternating order, as is the case 
for models without exchange and pair-hopping, then it is convenient to represent the configurations with non-zero trace by 
collections of segments. The weight of a segment configuration is determined by the length of the segments and the overlap 
between segments of different flavors (indicated by the hashed regions). 



In the simulation, one actually stores and manipulates M a , the inverse of the matrix defined in Eq. (|10|) . Fast 
matrix updates, similar to the ones detailed in Ref. ^3] allow to compute the new M a in a time 0(k^). The elements 
of this matrix also yield the measurement values for the Green function G a at the time intervals given by the operator 
positions (Tj for annihilation and rj for creation operators), 

G *( r ) = (g E M a (j,i)A(r,Ti -rj)}, (17) 

i,j=l 

A/ '\ J S(t-t') t' > . . 

A ^ T ^ = \-6(r-r'-f3) r' < ' ( 18 ) 

Angular brackets denote the Monte Carlo average. Other observables can be measured by computing a trace. For 
example, the mean particle number can be obtained as 

= / Tr[K loc {[3 - T 2n )Q 2n (T 2n ) . . . OMK^jT^h] \ 
11 \ Tr[K loc (f3-T 2n )0 2n (T 2n )...0 1 (T 1 )K loc (T 1 )} /' 1 ' 

where n is the number operator. 

A computationally expensive part of this procedure is the evaluation of the trace in the acceptance rate of Monte 
Carlo moves. In general, there are certain combinations of operators which always yield a zero trace and checking 
these conditions beforehand allows one to avoid unnecessary computations of the trace. 

Models which do not contain exchange or "pair-hopping" processes, so that -ffi oc and the ^-operators are diagonal 
in the flavor indices a, constitute a special case. For these models, the creation and annihilation operators for each 
flavor must occur in alternating order and as shown in Ref. 0] the "segment" representation, illustrated in Fig. [3 is 
an efficient way of specifiying all the configurations of non-zero trace. In this scheme, configurations are represented 
as collections of segments (one collection for each flavor) , whose start and end points coincide with the positions of 
the creation and annihilation operators. The weight of a configuration can be expressed in terms of the lengths of the 
segments and the overlaps between segments of different flavors. 

Care must be taken to prevent the system from being trapped in a state which breaks a symmetry of H\ oc when 
it should not be. For studying paramagnetic (para-orbital) phases, averaging the Green functions is sufficient. To 
study broken symmetry phases, the Green functions corresponding to different spin (orbital) states must be allowed 
to evolve independently and to obtain a symmetry unbroken state (e. g. above some critical temperature) it is then 
important that the Monte Carlo sampling explores the whole configuration space. To avoid un-physical trapping, 
we introduce "swap" -moves, which exchange the operators corresponding for example to up- and down-spins in a 
given orbital. Because the calculation of the new M a -matrices requires explicit matrix inversions, which are O(fc^), 
swap-moves are costly, but a relatively small number of attempts is enough to assure an ergodic sampling. 
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A. Self consistency and ordered phases 



The preceding sections described the solution of an impurity model specified by a local Hamiltonian and hybridiza- 
tion functions. In dynamical mean field theory, the hybridization functions are fixed by a self-consistency condition 
relating the impurity model Green function l|17|) to the implied lattice Green function. The precise form of the equa- 
tion depends on the specific dynamical mean field equation chosen, so a general equation cannot be written here. A 
crucial point is that the information concerning symmetry breaking is carried by the hybridization functions F and 
enters the problem via the self-consistency condition. H\ oc (and therefore the matrix forms of the creation and time 
evolution operators) retain their symmetry unbroken form. 

In this paper we use semicircular densities of states with (possibly orbital dependent) full bandwidths At a . The 
self-consistency condition for translationally-invariant states, including both paramagnetic states and states with 
ferromagnetic or ferro-orbital order is (the — r follows from the definition of F a in Eq. ijffil) 

F a (r)=tlG a (-r). (20) 

States with a broken translational invariance may also be studied. For example, for bipartite lattices with simple 
two-sublattice Neel order or (in the case of the models with two-fold orbital degeneracy) two-sublattice orbital order, 
the condition becomes 

F a (r)=tlG- a (-r), (21) 

where a denotes the opposite spin or the complementary orbital 

In subsequent sections we illustrate the formalism via study of two models in which exchange interactions play an 
important role: the Kondo lattice model and the Hubbard model with a two-fold orbital degeneracy. 



IV. APPLICATION I: KONDO LATTICE 



A. Overview 



In the Kondo lattice model, a local spin-1/2 degree of freedom, S, couples via a coupling constant J which may 
be either negative ( "ferromagnetic" ) or positive ( "antiferromagnetic" ) to electrons which reside in a single orbital, so 
that Eq. (j2J becomes 

H loc =-nYs + JS ■ l -^ a O a ^b- (22) 

a 

The Kondo impurity model, i.e. a single spin subject to a Hamiltonian -ffxondo + ^bath with -ffbath fixed (no self 
consistency equation) and characterized by a constant density of states p near the fermi level, has been extensively 
studied. The physics exhibits a profound dependence on the sign of the exchange constant J: for ferromagnetic J the 
coupling scales asymptotically to zero according to 

(J c ff(cj) ~ l/ln(w)) so that the asymptotic low temperature and low frequency behavior is that of free moments 
decoupled from the conduction electrons. On the other hand, for antiferromagnetic sign the problem scales to strong 
coupling, leading to the formation of a Kondo resonance and the dissolution of the spin into the bath of conduction 
electrons. 

Less is known about the lattice problem. We summarize here some results which are relevant to the half-filled 
case studied in this paper. For a classical spin the sign of J is irrelevant and for a bipartite lattice and particle-hole- 
symmetric dispersion the ground state is an antiferromagnetic insulator for all J |17| . The paramagnetic phase of the 
classical model is characterized by disordered spins, and may be an insulator at large J or a metal at small J. In the 
metallic phase the spin disorder implies a nonvanishing scattering rate at the fermi level. 

For S = 1/2 quantum spins, fewer results have been presented. It is generally believed that the half- filled, bipartite 
antiferromagnetic Kondo lattice exhibits a large- J Kondo insulator phase (the lattice version of the Kondo singlet 
behavior) whereas for smaller J a phase transition to an antiferromagnet occurs [lil Il9j . For the ferromagnetic 
side even less is known. A very recent study of the ferromagnetic Kondo lattice model at n ^ 1, based on the 
"equation of motion approach" which does not capture the Kondo scaling, reports a transition from a ferromagnetic 
to a paramagnetic state with increasing doping 
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B. Formalism 



We now turn to the specifics of the solution of this problem using the new method, H\ oc is diagonal in the basis 
of total particle number, total spin and z-component of total spin. If the particle number is or 2, then the spin 
state is just the state of the local moment, if the number is 1, the spin state is singlet (S) or triplet (T TOs! ) with given 
m z . The eigenstates may thus be labeled as shown in Tab.[IJ where the first entry is the number of electrons and the 



Eigenstates 


Ener; 


;y 


|i> = |o,T> 







!2> = |0,|) 







13} = \l,S) 






14} = |l,Ti) 


y- 


/' 


15} = |1,T ) 


y- 


i' 


16} = |l,T-i) 


y- 




17} = |2,T} 


-2fi 




18} = |2,I} 


-2n 





TABLE I: Eigenstates and eigenenergies for the local part of the Kondo lattice hamiltonian. The first entry labels the number 
of electrons and the second entry the spin state: either impurity spin j, J. if the number of electrons is or 2 or the total spin 
S (singlet) T m (triplet with m z = m) if n = 1. 



second entry refers to the spin state. The singlet state is defined as S = ^=(| t, I) — 1 1, t)), with the first entry the 
conduction electron and the second entry the local moment spin direction. In this basis, the time evolution operator 
is diagonal, K(r)\n) = exp(— E n T)\n), with eigenenergies E n listed in Tab.[3 The creation operators for spin up and 
down become the sparse matrices 
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(24) 



V 

With these operators, the sampling then proceeds as described in the previous sections. 

An important issue for simulations of interacting fermion problems is the sign of the different contributions to the 
partition sum. For the Hubbard model we noted the empirical absence of a sign problem in Ref. |l3j | . This absence 
of sign is not unexpected: the density-density interaction is essentially classical (no exchange) and other simulation 
methods do not give rise to a sign problem in this case. One might expect the situation in the Kondo lattice model to 
be worse, because it contains explicit exchange processes. Indeed as can be seen from Eq. i|24|) . the matrix elements 
for transitions into or out of singlet states can be negative. However, since these negative matrix elements always 
occur in pairs, the trace in Eq. l|15fl is not a source of sign problems. Negative determinants of the _F-matrices could 
in principle lead to negative weights (note that the exchange processes lead to operator orderings not found in the 
Hubbard model). Surprisingly, we do not find a sign problem in our simulations of the Kondo lattice either. For the 
parameters used in most of this investigation, (3t = 50 or 100, — 10 < J/t < 1 and densities per spin n < 0.98, the 
average sign is 1. Configurations with negative weight exist, but contribute negligibly little to the partition sum and 
are hence not generated. On some occasions, we measured average signs which differed from one in the sixth or seventh 
decimal place, but the converged solutions were usually not affected in any way by negative-weight contributions. 

A particularly attractive feature of the hybridization expansion approach is the fact that stronger interactions lead 
to a lower perturbation order, independent of the sign of J. In Fig. [3] we plot the distribution of perturbation orders 
p(k^) = p(fcj) corresponding to the converged solutions for different values of J/t and (3t = 50. While the distribution 
shifts in a way which is comparable to the one observed in the Hubbard model 01 for J < 0, the effect is even more 
pronounced for J > 0. For all parameter values considered in this study, the perturbation order remains reasonably 
low and thus allows an efficient Monte Carlo sampling. 
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FIG. 3: Distribution of perturbation orders p(k) = p(kf) — p(k\) for the ferromagnetic (left panel) and anti- ferromagnetic 
(right panel) kondo lattice models at half filling and inverse temperature fit — 50. Note the different J-ranges in the two 
panels. The mean perturbation order shifts lower as the coupling magnitude J is increased. For antiferromagnetic coupling, 
this effect is much more pronounced. 



C. Paramagnetic phase 

1. Overview and classical limit 

In this subsection we consider the behavior of the model in the paramagnetic phase (with magnetism suppressed by 
symmetrization of the Green function). For orientation, we first briefly discuss the physics of the classical core-spin 
model, in the paramagnetic phase. As noted above, in the classical model the sign of the exchange is irrelevant, and 
in the paramagnetic phase the spins are disordered and provide a static spin-dependent scattering potential for the 
electrons. In the dynamical mean field approximation to the classical spin model one finds that the self energy is 
E(u;) = J^q/Qq with 2 J e ff = J/2 the up-down energy splitting arising from the diagonal part of the exchange term 
in Eq. (|22|l and the mean field function Q^ 1 is given by 

Q-l 

GZ 1 (uj)=LU + fi-t 2 ql . (25) 

(Go 1 ) -J^ S 

At half filling (fj, = 0) and at the fermi level (ui = 0) this equation has two solutions: 

Go 1 = 0, (26) 

So" 1 = iyft J 2 cS - (27) 

Eq. l|27|l describes a metal (SroG ^ at the fermi level) with a self energy 

£ = ~i r~T~~Fy~ (28) 



2 

off 



which has a non-vanishing imaginary part, corresponding to scattering of electrons off the static spins. As J e g — > t 
the fermi level density of states vanishes and the scattering rate diverges. For \J c s\ > t the relevant solution is that of 
Eq. J2EJ), which is the uj — > limit of the expected insulating result Qq 1 ~ u> (describing an insulator with E ~ j-). 
In the rest of this section we present results for the quantum model, where the physics depends on the sign of J. 



2. Ferromagnetic J 

We begin with the ferromagnetic case. The left hand panel of Fig. 0] plots the converged Green functions for the 
ferromagnetic couplings J/t = —I, —3, —6, —8 and —10 at half filling and at the low temperature (3t = 50. It is 
apparent that for J/t > —6 the Green function is weakly dependent on J and exhibits the slow decay with time 
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FIG. 4: Left panel: local Green functions for the half filled ferromagnetic Kondo lattice model at the indicated exchange values 
and temperature. For J/t > —3 the computed Green functions are very close to the J = value and for J/t > —6 the long time 
behavior is characteristic of a metal. The exponential drop of the Green functions for J/t < —8 shows that the system becomes 
insulating at these large couplings. Right panel: dependence of the density on chemical potential. The smooth behavior for 
J/t > —6 shows the absence of a gap at half filling whereas a gap is clearly evident in the curve corresponding to J/t — —10. 
Some suggestion of a precursor to a gap is visible at J/t — —6. 
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FIG. 5: Imaginary time correlation function for the local moment calculated for ferromagnetic (left panel) and antiferromagnetic 
(right panel) couplings at half filling for the J values indicated. Solid lines show results for fit — 50. The dotted lines show 
results for f3t = 100 and J/t = —3 (ferromagnetic case) and J/t — 0.6 (antiferromagnetic case). 



characteristic of a metal. As the exchange coupling is increased, the system eventually undergoes a metal-insulator 
transition at a critical value between J/t = —6 and —8 (see also Fig. EJ), which is considerably larger in magnitude 
than the classical-model critical value J/t = —4 (J e g = — t). The right panel shows the dependence of the particle 
number per spin, n, on chemical potential for several J values; for small coupling, a smooth evolution is seen with no 
indication of a gap, whereas the opening of a gap is evident in the n(/^)-curve for J/t = —10. 

The left hand panel of Fig. shows the impurity-model spin-spin correlation function Css{ T ) — (Sz(0)S z (t)) 
calculated for the ferromagnetic Kondo lattice. An initial drop (with a J-dependent magnitude) is followed by a 
saturation to an almost temperature- independent value. In the classical model in the paramagnetic phase C(t) = 1/4 
independent of time. The saturation seen in the quantum calculation thus indicates that the long-time behavior of 
the spins is essentially classical, qualitatively consistent with the ferromagnetic Kondo scaling discussed in Eq. J23J. 
The combination of a paramagnetic state and a saturated (non-vanishing) spin-spin correlator implies the existence 
of annealed disorder in the spins, in other words the existence of zero frequency spin fluctuations. In particular, the 
saturation evident in the data for J/t = —8 and —10 shows that the charge gap seen in G(t) does not imply the 
opening of a spin gap. 

In the classical model, the spin disorder in the paramagnetic phase leads to a non-vanishing self energy at lu — > 
(either divergent, in the insulating phase, or finite, in the metallic phase). Fig.EJshows the self energies calculated for 
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FIG. 6: Left panel: imaginary part of the electron self energy E(itj„) for the half-filled ferromagnetic Kondo lattice model 
at fit = 50, J/t — —1, —3, —6 and —10. The metal-insulator transition which takes place between the last two values of J 
is obvious from the change in the low-frequency behavior. Right panel: expanded view of the low frequency behavior of the 
self energy in the smaller- J "metallic" phase. Dashed lines demonstrate an approximate power-law decrease of 9mS as ui — > 
with exponents 0.25 for J/t — — 1, 0.45 for J/t = —3 and 0.55 for J/t = —6. The solid line is proportional to the theoretically 
expected pH asymptotic behavior QttiE ~ l/(lnu n ) 2 . 



the ferromagnetic quantum model. For J/t = — 10, the system is insulating and S diverges, as in the classical case. 
However, for the smaller \J\ the self energy clearly vanishes as lu — + 0, a behavior quite different from that found in 
the classical case. 

The differences between the quantum spin-ferromagnetic coupling calculations and the results for the classical model 
have, we believe, a common origin, namely the decoupling of the carriers and spins at low energies (as found in the 
single-impurity model, Eq. 1)23(1 ). This is directly seen from the comparision of the spin-spin correlator (which shows 
classical spins) and the metallic phase self energy (whose vanishing at small frequency suggests no scattering at the 
fermi surface). This physics was already noted by Biermann and co-workers |2l| in a study of a related model. These 
authors argued that one could obtain the low frequency behavior of the electron self energy by combining the Kondo 
scaling Eq. lll'.'il) with the perturbative formula for the self energy to obtain T,(cu) ps pJ^ s (uj) ~ (1 + lnpa;) -2 . The right 
hand panel of Fig.[|J]shows an expanded view of the lower frequency regime of the metallic phase self energies. One sees 
that at the frequencies accessible to us the self energy is better fitted by a weak, J-dependent power law (the dashed 
lines correspond to the exponents 0.25 for J/t = —1, 0.45 for J/t = —3 and 0.55 for J/t = —6). In particular, except 
perhaps at J/t = — 1, the curvature of the numerical data is opposite to the curvature predicted by the one- impurity 
form. We suggest that the power law arises from an interplay between the one-impurity ferromagnetic Kondo scaling 
and the density of states renormalization due to J. In particular, for J near the critical value for the metal-insulator 
transition one expects a vanishing density of states. However, we note that the temperature range is insufficient to 
rule out a low-T crossover to the form proposed in Ref. |2l|. Further study of the frequency dependence of the self 
energy, and in particular a more precise characterization of the power law associated with the metal-insulator critical 
point, would be of great interest. 



3. Antiferromagnetic J 

The physics of the antiferromagnetic Kondo lattice is markedly different from that of either the classical spin or 
the ferromagnetic S = 1/2 Kondo lattice. The left panel of Fig. [7| shows the electron Green function calculated 
for several (small) J values and low temperatures. Comparison to Fig. Q] shows that for all J/t > 0.4, G(/3/2) falls 
below the value 4/(7r/3i) ss 0.0254 expected for a fermi liquid and approximately observed in the ferromagnetic case. 
Furthermore, as T is decreased G drops rapidly, suggesting the opening of the gap expected for a Kondo insulator 
[l8|. We believe that even the smallest J will eventually become insulating, but that the gap is too small to be seen 
on the temperature scales we have studied. The right hand panel shows that for J/t =1.0 and 1.2 a gap in the 
excitation spectrum is evident in the n(p) curve. Also, as expected in the presence of a charge gap, we find that the 
imaginary part of the self-energy diverges as w n — > (not shown) . 

The spin spin correlation functions for antiferromagnetic coupling are shown in the right hand panel of Fig. EI The 
correlations decay rapidly with time, consistent with the formation of a gapped Kondo insulating state. While the 



10 



<3 



0.01 



0.001 




0.65 




0.55 



FIG. 7: Left panel: thick lines show the local Green functions for the antiferromagnetic Kondo lattice model for J/t — 
0.4, 0.6, 0.8, 1 and inverse temperature /3t = 50. Thin dash-dotted lines correspond to fit = 100 and J/t = 0.8, 1. Right panel: 
density per spin plotted as a function of chemical potential. The data for J/t — 1.0 and 1.2 are consistent with the opening of 
a charge gap. 
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FIG. 8: Left panel: spin gap A s and charge gap A c extracted from fits of the function acosh(A(r — /3/2)) to the spin-spin 
correlation functions and Green functions obtained for /3t = 100 and the indicated values of J. Right panel: same data plotted 
as a function of t/J on a semi-log scale. The results are consistent with the expected small- J behavior In A ~ —1/pJ. 



exponential decay may not be obvious from the fit = 50 data, their fit = 100 counterparts (shown as an illustration 
for J/t = 0.6 by the dotted line) can be reasonably well fitted to a function of the form acosh(A s (r — (3/2)), From 
these fits we extract the spin gaps A s shown in Fig. El Also plotted are the charge gaps A c , which we obtained from 
an analogous fit to the Green functions. The variation of the gaps with J is very rapid and (as shown in the right 
hand panel of Fig. [HJ is roughly consistent with the theoretically expected behavior In A ~ —1/pJ at small J, crossing 
over to A ~ J for J > t. Remarkably, we find that the impurity model spin gap is less than twice the charge gap, 
with the ratio A s /A c decreasing through 1 as J is decreased. We understand this as a precursor of the magnetic 
state which would exist at small J and low T if magnetic order were not suppressed. However, we caution the reader 
that the spin gaps at the larger J-values are so large they are difficult to determine accurately, while the charge gap 
is uncertain at small J because the Green functions do not very nicely fit to a cosh-function. 



D. Magnetic Ordering 

We now show that our method correctly captures the magnetic ordering phenomena characteristic of the Kondo 
lattice. As in the previous sections, we specialize to half filling, bipartite lattices, and particle-hole symmetry. For 
orientation, we first review the known results for the classical-spin case. At half filling the classical model has 
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FIG. 9: Ferromagnetically coupled kondo lattice. Left panel: Green functions obtained for J/t — — 1 and fit = 10, 20, 30, 
40 and 50. A magnetic transition, setting in at fit ~ 30 is evident from the appearance of a difference between spin up and 
spin down Green functions. Right panel: staggered magnetization m — n-f — nj, as a function of temperature. Solid lines: 
S — 1/2 model, J/t = —1 and —2. Dashed lines: results from the classical model for J corresponding to the same diagonal 
spin splitting. 

antiferromagnetic order at all coupling strengths |l7j . At very small J, the classical transition temperature grows 
as ~ J 2 /t. It reaches a maximum around J/t ks 2 and for large J decreases as T$ ~ t 2 /J. In the quantum 
ferromagnetic case, we expect the Tjy ( J)-curve to retain essentially the same shape. In the quantum antiferromagnetic 
case we expect a quantum phase transition to a singlet phase for J larger than a critical value |l8l Il9| . 

We now turn to the results for the quantum model, beginning with ferromagnetic couplings. At half filling ferro- 
magnetism is never found to be stabilized, whereas the left panel of Fig.[3]sriows that with use of the antiferromagnetic 
self consistency condition a spin polarization (difference between up and down Green functions) becomes apparent for 
fit > 20 and J/t = — 1 . The spin polarization is associated with the formation of a gap, as can be seen from the rapid 
time decay of the lower-T Green functions (compared for example to the paramagnetic solution for fit = 50 in Fig.^J. 
Hence the ground state of the ferromagnetic Kondo lattice model is an antiferromagnctically ordered insulator. 

The right-hand panel of Fig.Elshows the staggered magnetization. Around T/t = 0.033, the staggered magnetization 
to = rir — n i for J/t = —1 starts to increase rapidly. We also plot data for a larger magnitude coupling J/t = —2, 
as well as results calculated for the classical model at couplings corresponding to the same effective spin splitting 
(dashed lines) j2^ . Surprisingly, in view of the ferromagnetic Kondo scaling, the critical temperatures in both models 
are comparable. While the magnetization onset in the quantum spin case is more rapid, the low-T saturation value 
is apparently lower. 

On the mean field level, one expects a continuous transition of the form m 2 ~ T c — T '. In the quantum case, we find 
m 2 (T)-curves which are roughly consistent with the linear behavior of the classical model, although the magnetization 
drops somewhat more rapidly near the critical temperature. The numerical data for J/t = — 1 might even hint a first 
order transition. In addition to the steep drop near T c , an essentially paramagnetic solution remains apparently stable 
for some range of temperatures below T c . However, a definite statement would require a more detailed investigation 
of the behavior near the critical point. 

In Fig. ^| we show the staggered magnetization of the antiferromagnetically coupled model as a function of J/t at 
several fixed temperatures. On the small J side a strong temperature dependence is evident, reflecting the strong J 
dependence of the Neel temperature at weak coupling. For J/t > 0.75 the fit — 40 data provide a good estimate of 
the T = result. At J/t > 1 the staggered magnetization rapidly drops to zero. This is the quantum phase transition 
to the singlet, Kondo insulator phase. We observe that this phase transition occurs at a J which is small relative 
to the bandwidth. The dashed line indicates the T = result for classical spins. In this case, no transition to a 
paramagnetic insulator occurs. The right hand panel again shows the magnetization as a function of temperature. 
For J/t = 1, magnetic order sets in around T/t — 0.077, which is noticeably higher than the transition temperature of 
the ferromagnetically coupled model or the model with classical spins. We attribute this to the growth in J implied 
by the antiferromagnetic Kondo scaling. On the other hand, due to the tendency to form singlets, the magnetization 
to = nf — ni for the antiferromagnetic model with quantum spins saturates at m rs 0.2, which is considerably smaller 
than the staggered magnetization of the corresponding ferromagnetic system. 
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FIG. 10: Staggered magnetization of the antiferromagnetically coupled Kondo lattice model (half filling, bipartite lattice, 
particle-hole symmetry). Left panel: staggered magnetization as a function of J/t for fit — 20, 40 and 80. There is an 
antiferromagnetic state at small coupling (for sufficiently low temperature) and around J/t = 1 a quantum phase transition to 
a paramagnetic insulator. The dashed line shows the T — result for classical spins. The right panel plots m — n-\ — n\ as a 
function of temperature. We find that the transition temperature is considerably higher than for ferromagnetic coupling and 
that the magnetization saturates at a smaller value. This smaller magnetization is the result of stronger quantum fluctuations 
(Kondo divergence) and singlet formation. 
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TABLE II: Eigenstates and eigenenergies for the local part of the 2-orbital model. The first entry corresponds to orbital 1 and 
the second entry to orbital 2. 



V. TWO ORBITAL MODEL 



For a second demonstration of the power of the method we consider here the two orbital model studied by other 
workers as a model for the orbital selective Mott transition. The local Hamiltonian is 

+ ^(C7'- J)n 1)(7 n 2iCT 

a — 1,2 cr a — 1,2 a a 

"• / Win,t^.^ 1 .t + ^2,1^2,1^1,1^14 + h - c -)- ( 29 ) 

We adopt the conventional choice of parameters, U' = f/ — 2 J, which follows from symmetry considerations for d- 
orbitals in free space and is also assumed to hold in solids. We consider semi-circular densities of states of bandwith 
4£i and 4^2 for orbitals 1 and 2, respectively, with a fixed ratio t%lt\ = 2, and furthermore restrict ourselves to the 
paramagnetic phase by averaging over spin in each orbital. 

The half- filling condition for this model is n — — I J and the 16 eigenstates and their energies are listed in 
Tab. [H] In this basis, the propagators K are diagonal, while the creation and annihilation operators for the different 
orbital and spin states become sparse 16 x 16 matrices. For a given spin, no more than two creation (or annihilation) 
operators may occur in a row and we check this condition before actually computing the trace. 
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FIG. 11: Converged Green functions for the 2-orbital model with semicircular density of states at half-filling. The ratio of band 
widths is £2/^1 = 2 and the temperature /3ti = 50. The exchange coupling is fixed as J = [7/4. For U/ti = 4 (i.e. U < C/f), 
both bands are metallic, whereas for U/ti = 8 (i.e. U > U2), both bands are insulating. For U/ti = 6, which lies in between 
Ui and U2, the narrow band is insulating, while the wide band is still metallic. Dotted lines show the non-interacting Green 
functions for fit = 50 and fit = 100. 




FIG. 12: Left panel: distribution of the orders p(k a ), a = 1, 2 for U/ti = 4 and U/ti = 6. The average order is lower for 
the narrow band and decreases with increasing interaction strength. Right panel: imaginary part of the self-energies for the 
metallic states in the narrow and wide band, showing the strong correlations in the narrow band for U/tj = 4 and in the wide 
band for U jt\ = 6. At U jt\ = 4, the wide band is only weakly correlated. 



An issue of debate in recent years has been the occurrence of an orbital selective Mott transition in such 2-orbital 
systems with Hund's coupling and different band widths, t\ 7^ £2- Using exact diagonalization [6j to solve the impurity 
problem, Koga et al. |23j found that the narrow band becomes insulating at a smaller coupling than the wide band. 
For semi-circular densities of states with a ratio of band widths t^/ii — 2, the critical couplings were found to be 
approximately U°/ti = 5.4 and E/J/ii = 7. On the other hand, in earlier work using QMC simulations, Liebsch [24| 
concluded that the transition takes place simultaneously in both bands. The QMC method should be more reliable 
than a ED calculation with a small number of bath sites, but the straight forward extension of the usual auxiliary 
field approach [5j suffers from a bad sign problem in the presence of spin flip and pair hopping processes, which were 
thus ignored in Ref. [24|. Arita and Held |lfj| have recently used a new type of Hubbard-Stratonovich decomposition 
0, which reduces the sign problem, and a projective QMC algorithm in their study of the two-orbital model. They 
found evidence for an orbital selective Mott transition in the presence of spin exchange, yet a single transition when 
merely the Ising component of the Hund's exchange was taken into account. Their estimate of t/f was consistent 
with the value obtained in Ref. |23j, while the projective QMC method did not allow to compute results at large 
enough couplings to estimate U^. Other recent works fliil 1251 l2fj | report the observation of two successive first order 
transitions and highlight the importance of taking the full Hund's coupling into account. It is therefore instructive to 
test our new algorithm on this example. 
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FIG. 13: Illustration of the relationship between an orbital selective state and the ferromagnetic kondo lattice model. The 
left figure shows the spin-spin correlation function in the narrow (insulating) band. The right figure shows the imaginary part 
of the self-energy for the wide (metallic) band. As in the ferromagnetic kondo lattice model, the correlation function in the 
insulating orbital saturates at large times, while the self-energy of the metallic band drops very rapidly as u) n — > 0. 

In Fig.^Jwe show converged Green functions for j3t\ = 50, J = U/4 and U/ti — 4, 6, and 8. The chemical potential 
corresponds to half-filling and we average over spin up and down in each orbital. Our continuous-time algorithm does 
not suffer from any sign problem at these parameter values and the hybridization expansion approach allows us to 
access large interaction strengths. Fig. ^] shows the distribution of orders in the two bands for the two coupling 
strengths U/ti = 4 and U /ii = 6. As expected, the average order is lower for the narrow band and the peaks shift to 
lower order as the interaction strength is increased. 

At U/ti =4, both bands are metallic, as can be seen in Fig. llll from the long-time behavior of G(r), which is close 
to the non-interacting solutions, shown by the dotted lines. The shift in the narrow band indicates that the electrons, 
while still itinerant, are becoming more strongly correlated. At U/t± = 8, both bands are insulating as follows from 
the rapid (approximately exponential) drop of G(r) to values much lower than the metallic solution. For U/t\ = 6, 
Fig. ^] shows that G(t) for the wide band is close to the non-interacting value, indicating metallic behavior, while 
G(r) for the narrow band drops exponentially, indicating an insulating state. 

The right hand panel of Fig. ^| shows the self energies of the metallic bands for U/ti — 4, 6. We see that for 
U/ti — 4 the wide band is weakly correlated (self energy small compared to frequency and to the bandwidth) but 
the narrow band is strongly correlated. At U/t% — 6 the wide band is strongly correlated (self energy larger than 
frequency and indeed comparable to half the bandwidth), but similar to a fermi liquid in the sense that ^sm'S(iuj) — > 
as to — > 0. However, a more detailed analysis reveals interesting differences with conventional fermi liquid behavior. 

As noted by Biermann et al. |2l|, the insulating orbital is effectively a local moment, which is coupled to the metallic 
orbital by the exchange coupling J. The usual Hund's rules imply that the exchange is typically of ferromagnetic 
sign; thus in the orbital selective phase one might expect the model to map onto a ferromagnetic Kondo-Hubbard 
lattice, with both an exchange coupling to a local moment and an on-site repulsion. Fig.ll3lshows that this is (at least 
qualitatively) indeed the case. The left hand panel plots the spin-spin correlation function of the insulating orbital. 
The initial drop and saturation behavior characteristic of the ferromagnetic Kondo lattice model is evident. (Note 
that because the orbital 1 can be empty or doubly occupied, the correlation function of the two orbital model at r = 
is slightly less than 0.25.) The magnitude of the initial drop is surprisingly large. The plot was made for U/t\ = 6 and 
J/U = 0.25 implying J jti — 0.75. In the Kondo model, J's of this magnitude lead to a much smaller decrease of Css 
from its initial value. For comparison, we also plot results for J/t = —0.8 in Fig. ^21 but as can be seen, an effective 
J/t w —6 is required to reproduce the 2-orbital results. Similarly, the right hand panel compares the calculated self 
energy of the wide band to the Kondo lattice self energy corresponding to a J/t chosen to approximately reproduce 
the drop in Css- The qualitative behavior with a rapid decrease at low frequency is the same in both models, but 
the quantitative agreement is not good, suggesting that much of the self energy of the metallic band arises from 
the U, rather than from spin-dependent scattering due to the Kondo coupling to the localized orbital. We have not 
yet run simulations at low enough temperatures to test the occurrence of the power-law behavior in the self-energy, 
demonstrated in Fig. for the Kondo lattice model. 

Our results indicate that the ferromagnetic Kondo-Hubbard model exhibits an interesting interplay between the 
on-site repulsive interaction and the Kondo coupling, leading to a much larger effective exchange coupling than implied 
by the bare parameters. Further exploration of this physics is an important open issue. 
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VI. CONCLUSIONS 

We have presented a formalism which extends a previously proposed diagrammatic QMC method to wide classes 
of impurity models. The idea is to expand the partition function in the impurity-bath hybridization function, while 
treating the local part of the Hamiltonian exactly. The resulting matrix formalism allows an efficient simulation 
of models with reasonably small Hilbert spaces. We have demonstrated the usefulness of the new approach with 
simulation results for the Kondo lattice and two orbital models. In both cases, the simulations in physically interesting 
parameter regions do not suffer from a sign problem. 

The new formalism opens up wide classes of questions for investigation. Systemtic investigations of quasiparticlc 
and magnetic properties of orbital selective Mott phases are now possible. We have provided direct calculations which 
support the conjecture of Bicrmann et al. 21] that the orbital selective Mott phase is in some qualitative sense 
described by an effective ferromagnetic Kondo lattice model, and we have further demonstrated that the Coulomb 
correlations in this phase play a very important role, leading to an effective coupling much larger than expected from 
the basic scales of the model. Concerning the Kondo lattice model, we have shown by comparing the ferromagentically 
and antiferromagnetically coupled cases that the renormalizations familiar from the one-impurity problem survive and 
have pronounced effects on the lattice problem, even at interaction scales of the order of unity. For example, the Neel 
temperature of the | J\ — 1 models differ considerably in the ferromagnetic and antiferromagnetic cases, which we 
believe is a result of the opposite renormalization of J in the two cases. For the ferromagnetic Kondo lattice model, 
we have discovered an unusual power law renormalization of the electron self energy which we propose is related to 
the density of states renormalization associated with the J-driven metal-insulator transition. Further investigation 
of this transition will be a fruitful subject for future research. For the antiferromagnetically coupled model we have 
located the Kondo-insulator to antiferromagnet transition and shown that the variation of the magnetization near the 
transition point is extremely rapid. 

Our method is from a conceptual and technical point of view appealing, because it does not require a double 
expansion in both the hybridization and the exchange couplings. The algorithm leads to manageable perturbation 
orders and, in the models studied so far, to an undetectably small sign problem in relevant regions of parameter 
space. In the presence of exchange processes, however, one has to compute the trace over all basis states in Eq. <|To|) 
explicitly. Because the Hilbert space grows exponentially with the number of orbitals, a straightforward application 
of the procedure introduced here becomes impractical for large impurity problems (a four site Hubbard cluster with 
256 basis states seems about the largest system one might want to consider). 

We see two possible ways to approach this problem. A straight-forward alternative is the above mentioned double 
expansion, which allows one to return to the economical segment picture |l3j | to represent the configurations and 
to devise efficient Monte Carlo moves which are compatible with the constraints of the model. Since the exchange 
couplings in many relevant models tend to be weak, the increase in the perturbation orders should still be manageable. 
What will happen to the sign problem remains to be seen. 

Another approach is based on the observation that most of the states in the exponentially large Hilbert space are 
of very high energy and are therefore not directly relevant to the physics. An important issue for future research is 
the development of "effective action" methods which will allow the elimination of high energy states, reducing the 
problem to one with a much smaller Hilbert space, to which the matrix formalism can be directly applied. 
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